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A fundamental assumption of the dynamical density functional theory (DDFT) of colloidal sys- 
tems is that a grand-canonical free energy functional may be employed to generate the thermody- 
namic driving forces. Using one-dimensional hard-rods as a model system we analyze the validity 
of this key assumption and show that unphysical self-interactions of the tagged particle density 
fields, arising from coupling to a particle reservoir, are responsible for the excessively fast relaxation 
predicted by the theory. Moreover, our findings suggest that even employing a canonical functional 
would not lead to an improvement for many-particle systems, if only the total density is considered. 
We present several possible schemes to suppress these effects by incorporating tagged densities. 
When applied to confined systems we demonstrate, using a simple example, that DDFT neccessar- 
ily leads to delocalized tagged particle density distributions, which do not respect the fundamental 
geometrical contraints apparent in Brownian dynamics simulation data. The implication of these 
results for possible applications of DDFT to treat the glass transition are discussed. 
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I. INTRODUCTION 

During the past decade of liquid-state research the 
classical dynamical density functional theory (DDFT) 
has proven to be a versatile and reliable tool for describ- 
ing the dynamics of interacting colloidal particles in a 
wide variety of situations. Building upon the success of 
equilibrium DFT [Tl [2], the dynamical theory enables 
first-principles calculation of the inhomogeneous density 
p{r,t) generated in response to a time-dependent exter- 
nal potential field T4xt(r, t) ?, . Within the DDFT frame- 
work a nonvanishing particle flux arises solely from gra- 
dients in a local chemical potential /z(r, t), derived from 
an equilibrium free energy functional. The theory has 
been successfully applied to many problems, including 
relaxation to equilibrium from a given nonequilibrium 
initial state [5] , the early and intermediate stages of spin- 
odal decomposition [3] and systems for which the time- 
dependence of K:xt(r,i) drives the system into nonequi- 
librium steady or stationary states [IHO] . In recent years 
the original formulation of DDFT has been extended to 
treat both systems and situations of ever increasing com- 
plexity. These more recent developments have lead to a 
better understanding of the influence of both nonpoten- 
tial fields (e.g. mechanical [TUHS], thermal [71II3]) and 
particle geometry (see e.g. [TS]) on diffusive colloidal dy- 
namics. 

Despite the undeniable success of the DDFT in ro- 
bustly capturing the qualitative features of p{r,t) for 
many problems of interest, the theory rests upon two fun- 
damental assumptions, both of which remain to be either 
systematically evaluated or improved upon. The first of 
these is the so-called adiabatic approximation, namely 
the assumption that the two-body correlations may be 
calculated from the instantaneous one-body density us- 
ing equilibrium statistical mechanical relations. The sec- 
ond major assumption is that the nonequilibrium chemi- 
cal potential /j,(r, t) generating the particle dynamics can 
be identified with the functional derivative of a grand- 



canonical free energy. Combining these two approxima- 
tions yields a closed equation for the one-body density, 
which does not require explicit knowledge of the higher- 
order correlations. 

In the present work we investigate the validity of em- 
ploying a grand-canonical density functional to treat 
many-body effects in DDFT. Problems can be antici- 
pated in confined systems with small particle number, 
for which the choice of ensemble strongly influences the 
equilibrium density profiles. More generally, use of the 
grand-canonical ensemble becomes questionable for sit- 
uations where the density field is strongly localized in 
space and contain only few particles. This can occur ei- 
ther as a direct consequence of minima in ycxt(r,i) or 
as a transient which may occur along the natural dif- 
fusive trajectory of p{r,t) through the space of density 
functions. Transient localization occurs quite naturally 
when considering the time-evolution of the density from 
sharply defined initial conditions, for which the positions 
of all A^-particles, or a subset thereof, are known pre- 
cisely. The individual density peaks of particles sharply 
located at i = remain well separated for short times 
and are normalized to unity (each peak contains a single 
particle). 

The above considerations become even more perti- 
nent when considering potential applications of grand- 
canonical DDFT to describe dynamic arrest and glass 
formation. By tagging the density field of a single particle 
in a dense hard-sphere liquid it has recently been shown 
[16] that the tagged density (the self part of the van Hove 
function) can exhibit a two step relaxation within DDFT, 
leading to dynamic arrest at sufficiently high volume frac- 
tions. Similar behaviour was also found for particles in- 
teracting via a Gaussian pair potential [17 . However, 
due to the use of an approximate free energy functional 
it remains unclear whether the observed slow dynamics 
arises purely from the presence of (possibly spurious) 
metastable minima in the free energy, or whether it is 
a genuine physical prediction of the DDFT which would 
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persist even when using an exact equilibrium functional. 
The work of Hopkins et al. |16[ I17j raises a fundamental 
question: Can a theory which involves only the one-body 
density field capture, even in principle, the transition to 
a nonergodic state? 

A major difficulty in answering the above question for 
realistic glass forming systems in two or three dimen- 
sions is the absence of an exact equilibrium free energy 
functional. It thus becomes difficult to disentangle errors 
in the dynamical theory from those induced by an ap- 
proximate free energy. For this reason we focus on the 
simplest nontrivial model for which the grand-canonical 
free energy functional is known exactly, namely a sys- 
tem of hard-rods in one dimension [THl HH] • Despite its 
simplicity, the hard-rod model has one feature in com- 
mon with glassy states occuring in higher dimensional 
systems: Both exhibit a partitioned phase space in which 
physically meaningful averages must be taken over a sub- 
set of phase space that is dynamically accessible. In 
glasses this partitioning occurs spontaneously as a ther- 
modynamic control parameter crosses some critical value, 
whereas for hard-rods the reduced phase space is always 
present, simply as a result of ordering the particles on a 
line. The confinement of a given rod by its neighbours to 
the left and right presents the prototypical glassy 'cage' 
and serves as a useful reality check for theories aiming to 
describe dynamical arrest in higher dimensions. 

The paper will be structured as follows: In Section 
\n\ we first will outline the microscopic dynamics under 
consideration, review the standard formulation of DDFT 
and define the functional to be employed. In Section 
|III| we consider the diffusion of a single rod and various 
methods by which grand-canonical contributions can be 
eliminated. In Section IV we identify the importance of 
tagging the individual density profiles. In Section |V] wc 
present numerical results for various test cases invlolving 
several interacting rods. Finally, in Section[Vl]we provide 
a discussion of our results and and identify some future 
challenges. 



II. FUNDAMENTALS 

A. Microscopic dynamics 

We consider a system of N colloidal particles in a time- 
dependent external potential V^xt(r, t) interacting via the 
spherically symmetric pair potential 4'{r). The total po- 
tential energy is given by 

i i<j 

As the individual colloidal trajectories are stochastic it 
is appropriate to adopt a probabilistic description of the 
particle motion. The probability distribution of particle 



positions is described by the Smoluchowski equation 



dPjt) 
dt 



0, 



(2) 



where P{t) = P({rJ,t). The fact that Eq.Q takes the 
form of a continuity equation expresses the conservation 
of particle number. Neglecting the influence of solvent 
induced hydrodynamic interactions (see [21| for a dis- 
cussion of the implications of this approximation) the 
probability flux of particle i is given by 



(3) 



with /3 = l/ksT. The force on particle i is generated 
from the total potential energy ([T]) according to = 
-V,C/jv. 



B. Dynamic Density Functional Theory 

Integration of Eq.([2]) over all but one of the particle 
coordinates leads to an exact coarse-grained expression 
for the one-body density profile 



dt 



-Vi-j(ri,i), 



(4) 



where the particle flux involves an integral over the 
nonequilibrium two-body density 

j(ri,i) =rfcBTVip(ri,i) + rp(ri,i)ViKxt(ri,t) 

+ ry"dr2Vi0(ri2)p(')(ri,r2,i), (5) 

with mobihty F = Do/fcsT. Eq.Q is the first in a hier- 
archy of equations for the n-body density. In equilibrium 
the flux is identically zero and Eq. Q reduces to the first 
member of the Yvon-Born-Green hierarchy |22| . 

The integral term in ([s]) may be approximated as an ex- 
plicit functional of the one-body density using the meth- 
ods of equilibrium DFT. The Helmholtz free energy is 
split into three contributions 

= J^id+ J^c^ + J dr I Vc^t{ri) piri), (6) 

where A is the thermal de Broglie wavelength and the 
ideal gas free energy is known exactly 

^id[p(ri)]= y'drip(ri)[ln(AV(ri))-l]. (7) 

The excess free energy functional Jox[p(r)] contains all 
information regarding the interparticle interactions and 
is connected to the average interaction force via the 
grand-canonical sum rule [T] 



J dr2 Vi(/i(ri2)/9(r2|ri) = Vi 



6p{ri) 



(8) 



where we have introduced the conditional distribution 
p('"2ki) = P^^H^1j''2)/p('"i), i-e. the average number 
density at r2 given a particle is fixed at ri. 
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The assumption that Eq.Q remains vaUd in nonequi- 
Ubrium is the so-called adiabatic approximation. It is 
equivalent to assuming that the pair density p^^^ (j^i , '"2 , 
relaxes instantaneously to the equilibrium pair-density 
corresponding to the current one-body density p{ri,t). 
The total particle flux ([5]) may thus be written in the 
form 

j(ri,i) - -rp(ri,t)ViM(ri,t), (9) 
where the nonequilibrium chemical potential is given by 



Sp{ri,t) 



(10) 



Combining Eqs.Q, ^ and ( |lO| ) yields the familiar form 
of the DDFT equation of motion 



dt 



Vi 



rp(ri,t)Vi 



Sp{ri,t) 



(11) 



In standard applications of DDFT the free energy T en- 
tering (11) is a grand-canonical quantity. Eq.(ll) thus 



describes the evolution of the density between two grand- 
canonical states, subject to the constraint that the aver- 
age particle number is conserved. 

For the present work the DDFT equation is solved nu- 
merically using finite difference techniques. The time in- 
tegration is performed using the Euler method and the 
spatial dervatives with central differences. The integra- 
tions required for the convolutions are evaluated using 
the trapezoidal rule, which in one dimension can be calcu- 
lated very efficiently due to the finite range of the weight 
functions. 



III. SINGLE PARTICLE DIFFUSION 



We begin by considering a single rod, whose diffusion 
is governed by the diffusion equation, resulting in a Gaus- 
sian density distribution. The diffusion equation can be 
recovered from the DDFT equation of motion ( [TT| ) in the 
case of vanishing interactions, i.e. J-cx = 0. In Figjljwe 
compare the results of DDFT (111 using the Percus func- 



tional with the exact Gaussian result for the relaxation 
of the density profile from the sharp initial condition 
p{x) = 6{x). The inset to Fig jl] shows the correspond- 
ing mean squared displacement (MSD) as a function of 
time. The DDFT reproduces neither the expected Gaus- 
sian form of the density profile nor the linear increase of 
the MSD as a function of time. Only at long times does 
the slope of the MSD approach unity as the local den- 
sity becomes very small for all x and the ideal gas term 
starts to dominate the free energy The deficiency 
of the theory lies in the nonvanishing contribution from 
the excess free energy, which leads to an effective self 
interaction of the density field and consequent enhanced 
relaxation rate. As originally pointed out by Marconi and 
Tarazona 3J, employing a grand canonical functional ef- 
fectively puts the system into contact with an unphysical 
particle reservoir, such that even for (N) — 1 the density 
distribution contains additional contributions from states 
with = 0, 2, 3, • • • . The interaction term thus does not 
vanish, as it should for a single particle, because states 
with N > I naturally incorporate interparticle interac- 
tions. The first step towards an improved theory is thus 
to eliminate, or at least reduce, the interaction between 
the physical particle and the reservoir. 



C. The Percus Hard-rod functional 

For our numerical investigations we will employ the 
exact excess Helmholtz free energy functional of hard- 
rods in one dimension [181 1191 . For an m component 
mixture of rods the functional is given by 

/oo 
dx (12) 
-00 

where the free energy density, $ = — no(x) hi(l — ni(x)), 
is a function of a set of weighted densities 

POO 

na{x) = J2 dx' p,{x')uj^^\x-x'), (13) 

where Pi{x) is the density profile of species i. The weight 
functions are characteristic of the particle geometry, with 
ujf\x) = 9(1 a;| - di/2) and u:f\x) = 5{\x\ - rfi/2)/2, 
where is the length of rod species i. 



A. Canonical correction 

In Refs. [23] and [Mj Gonzalez et al. proposed a scheme 
by which canonical equilibrium density profiles can be 
expanded in terms of grand-canonical averages. The 
method can be interpreted as a formal expansion of the 
canonical density profile in inverse powers of {N). Using 
this expansion, it was found possible to systematically 
correct the grand-canonical density profiles predicted by 
equilibrium DFT to achieve an improved description of 
a hard-sphere fluid conflned in a small spherical cav- 
ity. Generalizing the arguments presented in |23' '2T , the 
canonical average of an arbitrary function of the particle 
positions A{{ri]) may be expressed in terms of grand- 
canonical averages. To second order the expansion is 
given by 



{Af ^{A) + h{A) + f2{A) + -- - , (14) 
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FIG. 1. (Color online) The time evolution of density of a sin- 
gle particle for times t = 0.5, 1.0, 5.0, 10.0 (in units of / Dq). 
The length scale is set by the rod length d. We compare the 
exact result (grey dotted) with standard DDFT using the ex- 
act Percus functional (blue solid) and both first order (green 
short dashed) and second order (red long dashed) canonical 
corrections. The inset shows the corresponding mean squared 
displacement. 



where (•) is a grand-canonical average and (•)'^ is a canon- 
ical average. The correction terms are given by 

f,[A) ^ - \{{N ~ {N)fy^ 

h{A)^~\{{N-{N)f)f^ 



Due to the scaling of the partial derivatives in ( 14 1 the 
terms /i and /2 are of order (N)^^ and (N)^'^, respec- 

(ITl 



tively. The utility of the expansion ([14j) lies in its rapid 
convergence, even for very small values of (N) 23, 2A\ . 

Returning to the diffusion of a single particle, we now 
seek to use ( 14 ) to suppress unwanted grand-canonical 



contributions to the dynamics of p{r,t). At any given 
time we can use the instantaneous density profile p{r,t) 
to construct an effective external potential 

Kff(r, t) = Hpir, t)) - c(i) [p(r, t)] + ln(z), (15) 

where the fugacity term ln(z) is a physically irrelevant 
constant. The one-body direct correlation function 



5p{ri) 



(16) 



is evaluated using the instantaneous density. When 
employed in an equilibrium calculation, the potential 
VcG{r,t) will, by construction, yield the equilibrium re- 
sults for p(r, t) and all other quantities accessible to 
DFT. The adiabatic approximation is equivalent to as- 
suming that the higher-order nonequilibrium correlations 



are equal to the higher-order equilibrium correlations cal- 
culated in the presence of Vi2ff(r,t). 

We now define I{ri,t) as the average interaction force 
in the grand-canonical ensemble 



/(ri,t) = y"dr2Vi 



(17) 



where the subscript (gc) makes explicit the fact that the 
pair-density inside the integral is a grand-canonical aver- 
age. Using (|8| yields 



/(ri,t) = -fcBrp(ri,i)Vic«(ri,t) 



(18) 



Expressed in this way, it is rather natural to employ the 
expansion ( 14 ) to approximate the required integral in 



terms of known grand-canonical quantities 



dr2Wicf>{n2)p^^Hri,r2) ^ I{ri,t) + h{I{ri,t)) 



-h{I{rut)) 



(19) 



When calculated to all orders, the corrected interaction 
force should be zero for a single particle. We now pro- 
ceed to explicitly calculate the first two correction terms 
in ( |19[ ). The practical implementation of our scheme is 
as follows: At each time step in the numerical integration 
of ([2| the effective potential (151 is constructed from the 
instantaneous density. The partial derivatives required 
to evaluate the functions /„ appearing in ([l9]) are then 
calculated numerically by performing equilibrium DFT 
calculations in the presence of a fixed Vcsifiit). The se- 
ries ( 19 ) is then evaluated to the desired order and used 



to generate the density distribution at the next timestep. 
Grand-canonical contibutions to the time-evolution of 
the density p{r,t) can thus be suppressed on-the-fly in 
order to provide a more realistic description of the tra- 
jectory through the space of density functions. 

In Figjljwe show the time-evolution of the density of 
a single particle corrected using ( 19 1 to both first and 



second order. The series converges very rapidly and the 
second order results are virtually indistinguishable from 
the exact Gaussian function, despite the fact that (N) = 
1. Similar rapid convergence has been observed in the 
static case 1241. 



IV. MULTIPLE PARTICLE DIFFUSION 

We have now established that a systematic suppression 
of the grand-canonical contributions to the dynamics in- 
deed leads to improved results, at least for the trivial case 
of a single particle. However, application of the same 
procedure to systems with iV > 1 reveals an additional 
complication. 

Consider first the equilibrium average for a single infi- 
nite potential well in the canonical ensemble with = 1 
and in the grand-canonical ensemble with (A^) — 1 (see 
ngure ^ and b). The average value of a quantity is the 
average over the value of this quantity for all microstates 
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FIG. 2. (Color online) Schematic illustration of the mi- 
crostates that contribute to canonical and grand-canonical 
averages for particles in one and two infinite potential wells, 
(a) Canonical average A'' = 1 (b) Grand-canonical average 
(N) = 1 (c) Canonical average N — 2 (d) Grand-canonical 
average (TV) — 2. The second and third microstate in (c) are 
unphysical in the context of DDFT for localized particles and 
lead to a overestimation of the interactions. 



appropriate to the given ensemble. As a result, the av- 
erage density profile for canonical and grand-canonical 
case are different because of the additional microstates 
arising from coupling to an external particle reservoir. 

The situation is similar to DDFT applied to transiently 
localized particles virhich have had not had enough time 
to diffuse far avifay from their starting position. The av- 
erage interaction force is a grand canonical average, and 
the potential wells discussed above are manifest in the 
effective external potential (15). In the grand-canonical 



case the microstates with more than one particle in the 
well give rise to a nonvanishing interaction, which leads to 
the erroneous MSD. The canonical corrections suppress 
these unwanted microstates and the resulting treatment 
is closer to the canonical case. 

Next consider two infinite potential wells. The canon- 
ical correction series can again suppress the microstates 
with more than two particles in total. But the canonical 
ensemble still includes microstates with two particles in 
the same well. In the context of DDFT for two localized 
particles these microstates again lead to an unphysical 
interaction force, which does not vanish even if the two 
density peaks are far away from each other. The inter- 
action force on a particle resulting from number fluctu- 
ations in the Veff('",i) well generated by the particle can 
therefore be interpreted as a self interaction. 

For larger numbers of particles and wells each well is 
coupled to a reservoir formed by the other wells. As N 
increases the canonical and grand-canonical average be- 
come increasingly similar and the canonical correction 
decreases in magnitude, rendering it useless for the sup- 
pression of the self interaction. 
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FIG. 3. (Color online) (a) Test case of two overlapping Gaus- 
sian density peaks, (b) The grand-canonical average force 
Ji — — fcflTpi(ri)ViCi(ri) (green solid line) acting on the 
density field pi(ri) (right peak). Negative values indicate a 
force pointing to the left. Also shown are the full first order 
canonical correction to the average force (long dashed line) 
and the same quantity without the cross correction /i ^ (short 
dashed line), (c) The three terms contributing to the first or- 
der canonical correction. The length scale is given by the rod 
length. 



What is required is a method to prevent the averaging 
procedure from including microstates states with more 
than one particle per localized density peak. 



A. Tagged particle approach 

An improved description can be achieved by viewing 
the A^-particle system as an A^-component mixture, in 
which each component corresponds to a single particle. 
This allows to precisely locate particles, even after their 
densities started overlapping. 

The DDFT for an arbitrary m-component mixture is 
a set of m coupled equations 



dpi{ri,t) 
dt 



= Vi 



snip,}] 

Spi{ri,t) 



(20) 



where i labels the species. By tagging each individual 



particle Eq. ( 20 1 constitutes a set of m partial differential 
equations in D -f 1 dimensions (D space dimensions and 
time) for the time-dependent one-body density fields. On 
the other hand, the Smoluchowski equation ([2| for such 
a system is a partial differential equation in -|- 1 di- 
mensions for the probability distribution. A numerical 
solution of the Smoluchowski eqution is therefore very 
demanding and only possible for very few particles. 
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Secondly, in real applications one is not forced to tag 
each particle individually and ( 20 ) offers complete flexi- 



bility as to which subsets of particles are associated with 
distinct species. For example, in the dynamic test parti- 
cle calculations of [K] only a single particle was tagged 
and the rest of the system treated as a second component. 

The application of the full correction series with a in- 
finite number of terms to such a tagged system corre- 
sponds to a physically sensible average, that does neither 
suffer from particle fluctuations, nor from the combinato- 
rial effects discussed above. However, for large numbers 
of species and higher orders the canonical correction se- 
ries becomes rather unwieldy, but progress can be made 
by only keeping some of the terms. This can be demon- 
strated by considering a pair of interacting particles, each 
carrying a distinct species label. 

To first order, the canonical correction for such a bi- 
nary mixture is given by 

{Ar - {Ar + fi\A) + + f,\A), (21) 

where the correction terms are given by 



_ i djNo) d^{A) 
2 9/3^0 d(No)'^ 

d(No) d'^(A) 
dPfii d{No)d{Ni) ■ 



Higher order terms are of similar form, but involve higher 
derivatives. 

In analogy with the single particle calculation pre- 
sented in Section III the first order correction (21 1 could 



be used to determine the time evolution of two tagged 
density peaks. However, as argued in section |III A[ the 
correction of the interaction at each time is equivalent 
to an equilibrium situation. So we can assess the rela- 



tive magnitudes of the three terms in (21 1 by considering 



a static situation for which the two tagged profiles are 
fixed to be overlapping Gaussians (see FigjS^), although 
the detailed functional form chosen is not important. 

Choosing A = /i(r,t) from (18 1 as target quantity in 
Eq.(21|, where particle 1 is defined to be on the right. 



we numerically evaluate each of the three terms in the 
canonical correction series for a given separation between 
the Gaussian peaks (see FigjsJ:). The three terms can be 
interpreted as the corrections to the force acting on par- 
ticle one. fl^ corrects for grand-canonical fiuctuations 
in pi , for the increased interaction of po with pi due 
to fluctuations in po, and for the interaction between 
fiuctuations in pq and pi. If interactions are shortrange, 
it is clear that the term fl^ makes the dominant contibu- 
tion to the first order correction and that the small term 
fi^ may be neglected to a good level of approximation. 

Just keeping the terms and fl^ corresponds to 
a first order canonical correction for each component 
separately. As we have shown in section |III| this cor- 
responds to a first order suppression of the unphysical 
self-interaction, as each component consists of only one 
particle. Neglecting the mixed terms of the full canonical 



correction series leads to a system where each component 
does not interact with itself, but only with the other com- 
ponents. 

The conclusion which we draw from this simple test is 
that a suppression of the self-interaction in a system of 
N tagged density peaks, in which each peak represents a 
separate species and starts from sharp initial conditions, 
results in a reasonable approximation to the true Brow- 
nian dynamics at short and intermediate times. How- 
ever, the prohibitive numerical demands of performing 
the canonical transformation for more than a few parti- 
cles (expecially if higher orders of correction are required) 
make desireable an alternative scheme for eliminating the 
self interaction. 



B. Eliminating the self interaction 

1. Wtdom-Rowlmson model 

A well established model for which particles of the 
same species do not interact is the Widom-Rowlinson 
( WR) model [IS [55] . An approximate density functional 
for the m-component WR model has been derived [57] 
which provides a reasonable description of both the bulk 
phase behaviour and structural correlations. We propose 
to identify each species of an A^-component WR model 
with an individual particle. Employing the Widom- 
Rowlinson functional in this fashion guarantees the ab- 
sence of self interactions, but treats the interactions in 
an approximate fashion as the functional is not exact. In 
equilibrium with finite particle number (e.g. in confined 
geometries) the Percus and tagged WR functionals will 
yield different results. Due to the absence of unphys- 
ical self interactions the equilibrium density profiles of 
the tagged WR model should lie closer to the results of 
Brownian dynamics simulation than those obtained from 
the Percus functional. 

In one dimension the approximate excess free energy 
of the WR model is given by [27j 

/m 
dx^x) $ = 5]n^0. (22) 
1=0 

with weighted densities uq and for each component of 
the mixture. are the first derivatives of the OD excess 
free energy w.r.t. n\ and given by 

(^i = In ( 1 — m + cxp(zj) \ — Zi. (23) 

The fugacities Zi can be calculated from the implicit 
equation 

z, exp(zj) = ( 1 - m + ^ exp(zj) j nl (24) 
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FIG. 4. (Color online) Free expansion of five hard rods from 
a dense state at times t = 0.1,0.3,1.0. For clarity, curves 
for t = 0.1 and 0.3 have been shifted upwards by 1 and 0.5, 
respectively, and we show only positive x values. The rods 
are initially located at a; = 0,±2.5,±5, corresponding to a 
volume fraction around 0.4. Numerical results obtained from 



the Percus (12 1, Widom-Rowlinson (22 1 and subtraction (251 
functionals are compared with the results of Brownian dy- 
namics simulation. The simulation errorbars in this and all 
subsequent figures are smaller than the symbols themselves. 
The unit length is set by the rod length d, the times are in 
units of /Dq. 
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FIG. 5. (Color online) Evolution of the total density for the 
relaxation of four rods of length d — 1.6 in a periodic po- 
tential with minima at integer values and periodic boundary 
conditions. The lengthscale is given by the period h of the 
external potential, the timescale by h? /Dq. Shown are the 
results from simulations and DDFT calculations using the 
Percus functional ( 12 I and the subtraction functional ( 25 1 for 
, 100.0. 



times t = 1.0, 5.0, 25.0, 100.0. The corresponding evolution of 
a tagged density field is shown in Fig|6] 



The one particle direct correlation function required to 
calculate the average interaction force (181 is obtained 



by functional differentiation of the excess free energy. 



2. Subtraction of the self energy 

A simpler, albeit ad hoc method of suppressing the self 
interaction is to first calculate the excess free energy of 
full N component mixture using the Percus functional 
(12 1 and then subtract the individual excess free energy 
of each density peak. The remainder thus provides an 
approximation to the desired excess free energy arising 
from interaction between different species. This prescrip- 
tion amounts to employing the 'subtraction' functional 



(25) 



in the DDFT equation \2Q\. While the ansatz (25 1 is not 



justified from fundamental principles, it nevertheless has 
a certain physical appeal. In particular F^^^ vanishes for 
the case of a single particle and becomes exact for many 
particles in the low density limit. A similar approach was 
taken in |16j , in which an explicit self interaction term of 
a tagged particle was omitted from the Ramakrishnan- 
Yussouff functional [55] in order to recover the exact sin- 
gle particle diffusion. 



V. NUMERICAL RESULTS FOR SEVERAL 
RODS 

A. Free Expansion 

In order to compare the dynamics generated by the 
WR functional (22) with those generated by the sub- 



traction functional ( 25 ) we consider the free diffusion of 



five rods from sharp initial conditions. The initial delta 
function distributions are each separated by 2.5 particle 
diameters (corresponding to a volume fraction of around 
0.4). This choice ensures that the density fields have not 
entered the low density regime for the intermediate times 
at which overlap between neighbouring density peaks be- 
comes significant. The present test, which is similar to 
that considered inj3], thus constitutes a genuine test of 
the functionals ( [22] ) and (251 beyond the second virial 
level. In Fig|4] we show the time evolution of the total 
density 



(26) 



with m = 5, for three different times generated using the 
Percus ( 12 1, WR ( 22 ) and subtraction ( |25| functionals in 
the multi-component DDFT equation ( |20[ ). 

As was the case for an isolated particle(see Figjl]), the 
short time relaxation of the density peaks generated by 
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the Percus functional is too fast when compared with the 
Brownian dynamics simulation data. In contrast, the re- 
laxation of the total density predicted by both the WR 
and subtraction functionals is in almost perfect agree- 
ment with the simulation results. This good level of 
agreement supports our argument that the self interac- 
tion is the primary source of error induced by a grand- 
canonical generating functional, at least for short and 
intermediate times. Close inspection of the data reveals 
that the subtraction functional describes the simulation 
data slightly better than the WR functional, but the dif- 
ference is marginal. However, if the initial packing of 
the rods is denser, so that larger local densities occur, 
the WR functional becomes less reliable and predicts a 
unphysically fast relaxation, while the subtraction func- 
tional still gives a good description of the simulation data. 




-4-3-2-101234 



X 



B. Relaxation through a highly correlated state 

We now consider a further test case for which rods 
with a length of 1.6 relax from sharp initial conditions 
in a periodic external potential /3Kxt(a^) = — 2sin(27rx). 
The initial positions are chosen such that a particle is 
located at every second potential minimum. This situ- 
ation was originally suggested in [3] as a toy model for 
the study of slow dynamics. The external potential and 
rod length are constructed in such a way that two rods 
cannot be simultaneously at the bottom of neighbouring 
potential minima. Transport of particles between neigh- 
bouring minima thus requires a correlated motion of all 
A'' rods and it is the rarity of such events which leads 
to a long relaxation time. For long times an equilibrium 
solution is reached in which each potential well is equally 
populated. 

In FigjS] we compare the time evolution of the total 
density from standard DDFT employing the Percus func- 
tional ( [T2| ) and from a tagged particle calculation em- 
ploying the subtraction functional ( 25 1 with Brownian 
dynamics data. 

For short times the relaxation generated by the Percus 
functional is faster than that of the subtraction func- 
tional, in keeping with the intuition obtained from the 
case of a single free particle (see section III ). Surprisingly, 
for later times the relaxation of the subtraction functional 
profile overtakes that of the Percus profile, arriving more 
quickly to the equilibrium distribution in which, on aver- 
age, half a rod is located in each well. Use of the Percus 
functional thus yields better agreement with the simu- 
lation data than the, supposedly superior, subtraction 
functional. However, the good performance of the stan- 
dard theory arises from a fortuitous cancellation of er- 
rors, which can be appreciated by looking at the time 
evolution of a single tagged density peak. 

In Figjl we show the time evolution of the tagged den- 
sity corresonding to the third rod (labelling from left to 
right in Figjs]). The other tagged densities evolve identi- 
cally with time. Inspection of the profiles for t = 5 reveals 



FIG. 6. (Color online) A tagged density for same system 
considered in Fig[5] at t = 1.0,5.0,25.0,100.0. The Percus 
functional generates an unphysical growth of density peaks 
in next-nearest-neighbour minima which is suppressed when 
employing the subtraction functional. The inset shows the 
difference in the total density between neighbouring potential 
minima. Units as in Fig. |5] 



the strange behaviour of the tagged density from the Per- 
cus functional. Physically, it is reasonable to expect that 
the density peak will first diffuse into the neighbour- 
ing wells before spreading further to the next-nearest- 
neighbours, and so on. However, the Percus DDFT pre- 
dicts that the density in the next-nearest-neighbour well 
grows more rapidly than that in the neighbouring well. 
For intermediate times, the unphysically large density 
which has built up in the two next-nearest-neighbour 
wells then pushes back on the central peak and slows 
its decay (i.e. generates a component of the self interac- 
tion which tends to confine the remaining density in the 
central peak). 

In contrast to the Percus functional, the subtraction 
functional has a strongly reduced self interaction and 
does not suffer from an unphysical decay of the tagged 
density. Unfortunately, the improved description pro- 
vided by use of the subtraction functional destroys the er- 
ror cancelation presented by the Percus functional profile 
and thus relaxes much faster than the simulation data. 
The relative relaxation rate of the two approaches can be 
appreciated from the inset to Figj6] which shows the dif- 
ference in total density between neighbouring potential 
minima. 

The important message which emerges from this test 
case is that a qualitatively good description of the to- 
tal density relaxation is not sufficient to conclude that 
the DDFT is really capturing correctly the underlying 
dynamics of the tagged density profiles. 
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FIG. 7. (Color online) Equilibrium total density profiles of 
two hard rods confined to a slit with soft walls (Yukawa po- 
tential with screening length 0.2) located at positions x — ±5. 
The results are obtained from DFT Calculations with the Per- 
cus functional, the subtraction functional and the two com- 
ponent Widom-Rowlinson functional. The unit length is set 
by the rod length. 



FIG. 8. (Color online) The tagged particle density for the 
situation of figure [7| The exact density profile clearly shows 
how the particles are constrained to one side of the slit, as the 
particles cannot move through each other. At the wall, the 
exclusion zone for the left particle is clearly visible. The DFT 
profiles cannot capture this feature and the tagged particle 
densities of both particles are identical and delocalized over 
the whole slit. Units as in figure [7] 



C. Cage dynamics 

In two recent papers Hopkins et al. jl6[ I17j have 
shown that standard grand-canonical DDFT (in three- 
dimensions) can reproduce the two-step relaxation of the 
self part of the van Hove function [22] characteristic of 
glass forming systems. At sufficiently high volume frac- 
tions a divergent relaxation time was identified. In these 
studies the one-component Gaussian-core [T7] and hard- 
sphere [16' models were formally viewed as a two compo- 
nent mixture consisting of A^ — 1 particles of species d and 
a single particle of species s, thus tagging an arbitrary 
particle. A sharp initial condition was taken for particle 
s and the corresponding equilibrium distribution for the 
density field of the d-component 



Ps{r,0) 
Pd{r,0) 



Sir) 
Pb9{r), 



(27) 
(28) 



where pb is the bulk density. 

Interpretation of the results of [HI [17] was complicated 
by the fact that an approximate quadratic density func- 
tional was employed in the DDFT equation. It remains 
an open question whether the observed slow dynamics, 
which arise from a minimum in the equilibrium free en- 
ergy, are genuinely associated with the glass transition, 
or rather an indication of freezing within the theory. We 
think that some light can be shed on this issue by consid- 
ering the very simple case of two rods confined between 
impenetrable walls. The unique ordering of the rods on 
the line restricts the kinetically accessible phase space. 
As this is a characteristic feature of arrested states in 
general, we believe that a robust account of the equilib- 
rium tagged density distributions of two rods confined to 
a slit is a neccessary pre-requisite for any DDFT claim- 



ing a realistic description of dynamic arrest in higher- 
dimensional systems. 

In FigjTjwe compare the equilibrium total density from 
DDFT with simulation data for two rods confined by the 
potential 



a;— 4.5 



4.5 — an 



/3Kxt(a;) 



4.5 4.5 



(29) 



with decay length I = 0.2. This gives a volume fraction 
around 0.2. For this situation it can easily be shown 
that the canonical equilibrium density of the left particle 
is given by 



x+d 



(30) 



The primary observation to be made from FigjTjis that 
the Percus functional generates a packing structure at the 
wall with a well developed peak and a dip. The dip is ab- 
sent in the exact profile and the peak is less pronounced. 
It should be recalled that the Percus functional profiles 
represent the exact grand-canonical result for a slit with 
(TV) = 2. The peak is slightly underestimated by the WR 
functional ( 22 ) and entirely absent from the subtraction 



functional profile. Overall the WR functional shows the 
best level of agreement with the simulation data for the 
total density p{x). 

On the basis of the total density profiles shown in FigjT] 
it is tempting to conclude that the WR functional pro- 
vides an acceptable description of the equilibrium den- 
sity distribution. Indeed, this is true if one is inter- 
ested only in the total density p{x). A completely dif- 
ferent picture emerges, however, when considering the 
two tagged density profiles. In FigjSjwe show the tagged 
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densities corresponding to the same situation as shown 
in FigjTj The difference between the theoretical predic- 
tions and the simulation data is dramatic. All of the 
DFT approaches yield identical equilibrium tagged pro- 
files pi{x) = P2{x) — p{x)/2, clearly demonstrating that 
the DDFT time evolution does not respect the fact that 
hard-core particles cannot pass through each other. At 
some point during the time evolution the density fields 
always tunnel through each other to arrive at unphysical 
regions of phase space. Given that none of the functionals 
considered in the present work are capable of capturing 
this most elementary 'caging' dynamics, we find it un- 
likely that the DDFT employed in [HI [TT] is capable of 
capturing a true nonergodic transition. The constrained 
order of the rods is encoded in the hierarchy of correla- 



coUective motion of rods (see section V B ) demonstrated 



tion functions (see Section Il||II B ) in a complicated way, 
so one can not expect that truncating this hierarchy pre- 
serves this property. 



VI. DISCUSSION 

We have analysed the shortcomings of DDFT for the 
description of localised density distributions using a sim- 
ple one-dimensional model of hard-rods. This model has 
the advantage that the exact grand-canonical equilib- 
rium functional is known, thus removing the possibility 
of dynamical artifacts arising from use of an approximate 
functional. It is well known among DDFT practitioners 
that the standard theory predicts a relaxational dynam- 
ics which is systematically too fast when benchmarked 
against Brownian dynamics simulation. By employing a 
formal expansion of the canonical average we conclude 
that this deficiency is due to the unphysical self interac- 
tion arising from grand-canonical coupling to a reservoir. 

In situations where few particles are strongly confined, 
the presence of a self interaction leads to equilibrium den- 
sity profiles which differ from thise obtained in Brownian 
dyanmics simulation. We anticipate that the self inter- 
action will become relevant for large systems in cases for 
which large unbalanced forces are present, e.g. during 
the relaxation of large density gradients in a nonequilib- 
rium profile. The corrections in the force resulting from 
removal of the self interaction will modify the relaxation 
timescale and slow the relaxation relative to standard 
grand-canonical DDFT. 

The only method by which the self interaction can be 
removed is to individually tag the density field of each 
particle and then employ either the canonical expansion 
series, a Widom-Rowlinson-type model or ad hoc sub- 
traction of undesirable contributions to the excess free 
energy functional. Of these three possibilities, the lat- 
ter proved to be both the most reliable and simplest to 
implement. 

Once the self interaction has been dealt with appro- 
priately, we find that the free expansion of any number 
of interacting rods can be well described using DDFT. 
However, the considered test case which focused on the 



that the previously reported good agreement between 
DDFT and simulation [3] is due to cancellation of errors. 
Removal of the self interaction corrects the unphysical 
aspects of the tagged density relaxation, with the con- 
sequence that the results for the relaxation of the total 
density become worse (even faster than standard DDFT) . 
The fact that DDFT clearly has difficulty in describing 
the relaxation through highly correlated states raises sus- 
picions about its ability to capture the slow dynamics of 
dense systems. 

Although our one-dimensional model provides only a 
crude description of a real fiuid, it nevertheless enables 
some of the subtleties associated with a partitioned phase 
space to be investigated within a simple setting. Moti- 
vated by the recent work of Hopkins et al. [ini E] we 
presented the case of two confined rods as a toy model for 
the cage effect, whereby particles in a glass are localised 
by the geometrical contraints of their neighbours. Given 
that the very simple geometrical contraints on the tagged 
particle densities are not respected (see Figjs]) it would 
be remarkable if the same theory turned out to be capa- 
ble of describing the spontaneous localization associated 
with the glass transition in two- and three-dimensions. 
One possibility is that the particular functional employed 
in (161 117] , namely the quadratic Ramakrishnan-Yussouff 
functional [SH], is somehow able to compensate for the 
errors arising from a grand-canonical treatment of the 
density. However, such a scenario would seem to require 
a highly nontrivial cancellation of errors. 

By considering test cases in which finite numbers of 
rods are confined to a slit we now strongly suspect that 
any theoretical approach which is closed on the level of 
the one-body density, such as DDFT, will be unable to 
describe localization of the tagged density fields (at least 
when employing the exact equilibrium functional) . When 
working solely with the one-body density, effective inter- 
actions between the average quantities Pi{v) are imple- 
mented, whereas Brownian dynamics considers interac- 
tions between the density operators /5i(r) = i5(r — r^) 
before taking the average. This mean-field treatment of 
the interaction between tagged density fields appears to 
make the tunnelling of tagged density into geometrically 
forbidden regions unavoidable. We note that these lim- 
itations do not neccessarily pose a problem for the es- 
tablished density functional theory of freezing P^^. In 
such calculations the crystal is identified as a periodic 
variation of the total density field and no claim is made 
to identify any specific particle: the tagged densities are 
fully delocalized throughout the sample. 

To arrive at a theory which respects the geometrical 
constraints on the tagged density fields it is neccessary 
to go beyond a density based description and consider 
explicitly the dynamics of the higher-order density corre- 
lators, i.e. going beyond the simplest adiabatic approx- 
imation. Indeed, the importance of improving upon the 
standard adiabatic approach was recently identified by 
Haataja et al. [30 as one of the most important open 
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problems in DFT. A fundamental question is whether 
a finite order truncation of the dynamical hierarchy (of 
which ([4]) and (|5| constitute the first member) can cap- 
ture tagged density localization. 

Integration of the Smoluchowski equation over all but 
two of the particle coordinates leads to an equation of 
motion for the pair density involving a weighted integral 
over the three particle correlations. Specifically, one is 
required to evaluate integrals of the form 



p(2)(ri,r2) 



(-Vi0(ri3)). 



(31) 



The factor in square brackets can be identified as the 
conditional probability to find a particle at r^, given that 
there are particles at positions ri and r2- The integral 
thus represents the average force acting on a particle at 
ri due to interactions mediated by the other particles (an 
analogous integral provides the indirect force on a parti- 
cle a t r2 ). Making an adiabatic approximation, the inte- 



gral (|3l[) may be replaced by Vic^rliri), where c'^r^iri) 



is the one-body direct correlation function at ri calcu- 
lated in the presence of both the physical external field 
and a test particle fixed at r2 (hence the parametric de- 
pendence upon r2). This leads to 



9/j(^^(ri,r2) 
Tdt 



Vi 
V2 



f 5F 



\^P{ri)/2 



\5p{r2))i 



.(32) 



Combining Eqs.Q, ([5| and (32 1 leads to a closed theory 
for the dynamics of the one- and two-body density in the 
form of a pair of coupled first order (in time) differen- 
tial equations. A similar equation was employed in [3T] 
(see their Eq.(9)) in which the three-body density was 
treated using a superposition approximation. Crucially, 
in ( 32 1 the free energy functional to be differentiated con- 



tains an external field consisting of the physical external 
potential of interest, Vi3xt(''"), plus a 'test' particle held 
fixed at either position ri or r2, as indicated by the sub- 



script on the functional derivative. Equation (32 1 is still 



'adiabatic', in the sense that three- and higher-body cor- 
relations are determined from the nonequilibrium p(r,t) 
and p'^2)(ri, r2, i) using equilibrium statistical mechan- 
ical relations. Nevertheless, this extended theory goes 
beyond (11), as the pair density is no longer tied to the 



instantaneous value of the density and relaxes instead 
on a finite timescale. Note that the formal elimination of 
p(2) (n, r2) from the coupled equations generates, in prin- 
ciple, an equation for p{ri) alone which includes memory 
effects j32|- In practice, however, the nonlinearity of the 
equations does not allow an explicit form for the memory 
kernel to be obtained and memory effects remain implicit 
to the coupled system of equations. 

In equilibrium, Eq.(32| predicts that the pair corre- 



lations are those obtained from a test-particle calcu- 
lation using the chosen free energy functional. The 
exact single particle dynamics are recovered using the 
initial condition p*-^^(t'i, r2, 0) = 0. More generally, 
the correct normalisation of the initial pair density, 
J drij dr2 p^^''{ri,r2,0) = N{N — 1), is preserved 
throughout the time evolution, which is not the case 
when applying the standard adiabatic approximation, i.e. 
calculating equilibrium pair correlations in the presence 



of the effective potential ( 15 ). Although we defer a more 
detailed investigation of ( 32 i to future work, preliminary 



calculations indicate the results for the nonequilibrium 
one- and two-body density of hard-rods are considerably 
improved, relative to standard DDFT. Of particular in- 



terest will be the predictions of ( 32 1 for the density profile 



of confined systems in the limit of long times. 

Apart from our preliminary investigations of ( 32 ) the 
only explicit test of the adiabatic approximation of 
which we are aware was performed using continuous- 
time Monte-Carlo simulations of the Potts model sub- 
ject to Glauber dynamics [33j . In this study, initial sim- 
ulation configurations were prepared which reproduced 
the equilibrium one-body occupation number profile, but 
with nonequilibrium correlations between the occupation 
numbers. During the simulated time-evolution the relax- 
ation of higher-order correlation functions caused the the 
one-body profile to drift first out of equilibrium, passing 
through a sequence of transient intermediate states, be- 
fore returning back to equilibrium at long times. These 
findings clearly cannot be reproduced by theories based 
on p(ri,t) alone, as no distinction can be made be- 
tween states with the same one-body profile but different 
higher-body correlations. This issue may prove to be im- 
portant when considering systems with slow dynamics 
for which the one-body density remains constant during 
gradual structural relaxation processes. 

We thank M. Schmidt and A.J. Archer for stimulating 
discussions. Funding was provided by the Swiss National 
Science Foundation. 
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